% calculating r0 for given turbulence layer
% references in <<Principles of adaptive optics>>
function r0 = r0_generate(Cn2,lambda,H_bot,H_top,sp,L)
% lambda : wavelength [m]
% Cn2 : [km]
% H_top/bot : Computing distance with dimension [km]
% r0: The OUTPUT with standarized dimension [m], but not [cm]
% plane wave
% 
% if max(H_top) < 99 % dimension transforms from [km] to [m]
%     H_top = H_top * 10^(3);
%     H_bot = H_bot * 10^(3);
% end

k = 2*pi/lambda;

% 默认平面波
if nargin == 4
    sp = 'p';
end

if sp == 'p'
    if (strcmp(class(Cn2),'function_handle')) && (length(H_top) == 1)
        fun = @(z) Cn2(z*10^(-3));
        fun_itgl = integral(fun,H_bot,H_top);
        r0 = (0.423 * k^2 * fun_itgl).^(-3/5);
    else % when the inputted Cn2 is a Constant
        delta_z = H_top - H_bot;
        r0 = (0.423 * k^2 * Cn2 * delta_z).^(-3/5);
    end
end

if sp == 's'
    if (strcmp(class(Cn2),'function_handle')) && (length(H_top) == 1)
        fun = @(z) Cn2(z*10^(-3)).*(1-z/L).^(5/3);
        fun_itgl = integral(fun,H_bot,H_top);
        r0 = (0.423 * k^2 * fun_itgl).^(-3/5);
    else % when the inputted Cn2 is a Constant
        fun = @(z) (1-z/L).^(5/3);
        fun_itgl = integral(fun,H_bot,H_top);
        r0 = (0.423 * k^2 * Cn2 * fun_itgl).^(-3/5);
    end
end

% end of the function
end




